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Abstract 



We analyze the asymptotic behavior of a partial differential equation (PDE) model 
for hematopoiesis. This PDE model is derived from the original agent-based model 
formulated by Roeder et al. in [35J, and it describes the progression of blood cell 
development from the stem cell to the terminally differentiated state. 

To conduct our analysis, we start with the PDE model of |20j, which coincides 
very well with the simulation results obtained by Roeder et al. We simplify the PDE 
model to make it amenable to analysis and justify our approximations using numerical 
simulations. An analysis of the simplified PDE model proves to exhibit very similar 
properties to those of the original agent-based model, even if for slightly different 
parameters. Hence, the simplified model is of value in understanding the dynamics of 
hematopoiesis and of chronic myelogenous leukemia, and it presents the advantage of 
having fewer parameters, which makes comparison with both experimental data and 
alternative models much easier. 



Key-words Age-structured equations, hematopoiesis, chronic myelogenous leukemia, 
model simplification. 

1 INRIA Paris-Rocquencourt, BANG, BP105, F78153 LeChesnay Cedex. Email: marie. doumic- 
jauffret@inria.fr 

2 Laboratoire J.-L. Lions, Universite Pierre et Marie Curie, CNRS UMR7598, BC187, 4, place Jussieu, 
F75252 Paris cedex 05. Email: benoit.perthame@upmc.fr 

department of Mathematics, University of Utah, Salt Lake City, UT 84112-0090, USA. Email: 
kim@math.utah.edu 

* Contributed comparably to this paper. 



1 



1 Introduction 



Chronic myelogenous leukemia (CML) is a cancer of the blood and bone marrow that results 
in the uncontrolled growth of myeloid blood cells. More than 90% of all CML cases are 
associated with a gene abnormality, known as the Philadelphia (Ph) chromosome [36J. In 
addition, CML is highly responsive to treatment by the drug imatinib that specifically targets 
the gene abnormality [T7) . 

Recently, CML has been the focus of several mathematical models, which we summarize 
below. These works were motivated by the desire to explore the mechanisms that control 
the disease with the hope that this will lead to new therapeutical strategies. In 1991, Fokas 
et al. formulated the first mathematical model of CML [19]. In 2002, Neiman presented 
a model that accounts for the immune response to CML |30j. This work attempted to 
explain the transition of leukemia from the stable chronic phase to the accelerated and acute 
phases. A more recent work of Moore and Li in 2004 aimed to identify the parameters that 
control cancer remission [29]. Their main conclusion was that lower growth rates lead to a 
greater chance of cancer elimination. In 2005, Komarova et al. used methods of stochastic 
networks to study drug resistance with a particular view toward imatinib [23]. In the same 
year, DeConde et al. proposed a model for the interaction between the immune system and 
cancer cells after a stem cell (or a bone- marrow) transplant [13]. The main result of that 
work was that a slightly elevated autologous (pretransplant) immune response greatly favors 
remission. 

A new paradigm of cancer development emerged from the idea of cancer stem cells [5]. 
This hypothesis states that a variety of cancers originate from a self-replenishing, cancer 
population, now known as cancer stem cells. Using this idea, Roeder et al. developed a 
mathematical model of CML stem cells [35] . In their model, leukemia stem cells continually 
circulate between proliferating and quiescent states. This formulation contrasts with the 
alternative paradigm of Michor et al. in which leukemia cells differentiate progressively from 
stem cells to differentiated cells without circulating back to previous and more dormant 
states [28l [27]. Both the Roeder and Michor models are directed to studying the dynamics of 
imatinib treatment. However, each model also presents a general paradigm for hematopoiesis 
that describes blood cell development with or without CML. In this paper, we focus on the 
model of Roeder et al. 

The hypothesis of the Roeder model is based on experimental studies that demonstrated 
that within a two week period, nearly all hematopoietic stem cells enter cell cycle [HI ITT] . 
A summary of the experimental results and interpretations is also presented in [33]. Based 
on these findings, Roeder et al. formulated an agent-based model (ABM) to capture prob- 
abilistic effects and intrinsic heterogeneity of the stem cell population [35]. However, since 
the ABM is computationally demanding, Roeder et al. developed an analogous PDE model 
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that is based on their original ABM model [31] . In a parallel and independent work, Kim et 
al. also developed an analogous PDE model [23] . The main difference between the two PDE 
models, is that the version of Roeder et al. simplifies the original formulation by eliminating 
the explicit representation of the cell cycle, whereas the version of Kim et al. includes the 
full complexity of the original ABM. 

In this paper, we take the PDE model in [23] and simplify it as much as possible without 
altering the fundamental assumptions of Roeder et al. [35] . Then we conduct an analysis 
of the asymptotic behavior of the simplified model for hematopoiesis. The assumptions of 
Roeder et al. include the fact that stem cells exist in two growth environments: proliferating 
and quiescent. In addition, proliferating stem cells gradually progress toward further levels 
of differentiation until they differentiate completely. However, proliferating stem cells can 
reenter the quiescent state, in which they cease dividing and regress toward more primitive 
levels of differentiation, even up to the fully immature state. One can interpret that while 
quiescent stem cells cease dividing, they begin performing other regenerative functions such 
as returning to the stem cell niche as suggested in [35] . 

The paper is organized as follows. In Section (2J we present the PDE model from [23J, 
which proved to coincide very well with the simulation results of Roeder et al. in [3J5J- In 
Section [31 we present four approximations of the PDE model, one of which is very similar to 
that presented in [31] . In Section HJ we compare numerical solutions of the four approxima- 
tions with simulations of the original ABM. In Section [5j we analyze the asymptotic behavior 
of the simplest models given by Approximations 3 and 4. 



2 PDE version of Roeder's model 



In this section, we present a PDE version of the original ABM for hematopoiesis developed 
by Roeder et al. in |35j. This PDE model is taken from [23]. A state diagram for the model 
is shown in Figure I2TT1 

In this model, hematopoietic stem cells (HSCs) operate in two growth compartments: 
quiescent (Alpha) and proliferating (Omega). Stem cells transfer from Alpha to Omega or 
from Omega to Alpha at rates u and a, respectively, where 

a(x,A) = e-^f a (A)^ 
u(x,Q) = a wiB .e*' IB f w (Q). 

The coordinate x is a state variable that ranges from to 1, and the parameter a m i n = 0.002 
as estimated in [35J. The variables A and Q denote the population of cells in the Alpha and 
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01 Pr^: = -u(x^)A(x,t)+a(x,A)[ I n(x,c,t)dc+l Gl (x)Q*(x,t)) (2.2) 



Cl 



Omega compartments, respectively, and f a and f u are sigmoidal functions, whose definition 
can be found in Appendix [A1 

The equations for the PDE model in [23] are as follows: 

dA dA ' ^2 

~dt 

dA* — 

— = p r A(0,t) -u(0,tt)A*(t) (2.3) 

dn on on . 

~&t + Pd ~dx~ + lte = ~ a ^ Xl A )W 2 )( C M^ c ' *) ( 2 - 4 ) 
dfl* dQ* — 

-gf + Pd-fa = -a{x,A)l Gl {x)tt*{x,t) (2.5) 
with boundary conditions 

A(l,t) = 0, (2.6) 
Q(x, 0, t) = Q(x, c 2 , t) + lu(x, U)A(x, t), (2.7) 

n*(0, t ) = ^^-A*(t), (2.8) 
Pd 

Q(x,cf,t)=2Q{x,c];,t), (2.9) 

n*((kp d c 2 + c 1 ) + ,t)=2Sr((kp d c 2 + c 1 )-,t), k = 0,1,2,.... (2.10) 

In these equations, x G [0, 1], c G [0, C2], and 

Gi = ( U [^ dC2 + Cl ' ( fc + ) fit ' 1 1- 

Wo / 

The parameter c\ represents the duration of the combined S/G2/M-phases, and C2 represents 
the duration of the entire cell cycle. The set Gi is the set of x-values for which Q* cells are in 
the Gi phase (defined by c G [ci,c 2 )). In addition, the total Alpha and Omega populations 
are given by 

A(t) = [ A(x,t)dx + A*(t), 

(2-11) 







fi(t) = / / VL(x } c } t)dcdx+ I tt*(x,t)dx. 
Jo Jo Jo 



As shown in the state diagram in Figure l2"7Tj the Alpha cells are composed of two subpopu- 
lations: A and A*. As time progresses, A cells decrease their x-coordinate at rate p r , until 
they attain the minimum x- value of 0, at which point they enter A*. 

In a similar manner, the Omega cells are divided into two subpopulations: Q and Q*. 
The subpopulation Q corresponds to cells that have entered the Omega state from A, and 
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the subpopulation Q* corresponds to cells that have entered from A*. As time progresses, all 
Omega cells increase their x-coordinate at rate pa, until they attain the maximum x-value of 
1, at which point they differentiate into proliferating precursors. In addition, all Omega cells 
progress through the cell cycle as they undergo proliferation. The Gi phase corresponds to 
the period during which the cell generates new organelles. Only cells in the Gi phase can 
transfer to A. The other phases, S, G2, and M, correspond to the period during which a cell 
is undergoing mitosis. Cells in these phases cannot transfer. A cell's position in the cycle is 
measured by its c-coordinate, which ranges from to c-i. 




X = A Pr X = 1 



(sR) Q(x, ) and Q*(x) cells in Gi-phase transfer to A(x) 

Figure 2.1: State space for the PDE model of CML dynamics in [23]. The variable A(x,t) 
represents stem cells in the Alpha (nonproliferating) compartment with intermediate re- 
values, and the variable A*(t) represents cells in the Alpha compartment that have attained 
the minimum x- value of 0. The variables Q(x,c,t) and Q*(x,t) represent the stem cells in 
the Omega (proliferating) compartment, where Q and Q* correspond to cells transitioning 
from the A and A* subpopulations, respectively. 

The first term on the RHS of (12. 2p accounts for the cells that transfer out of A into Q. 
The transition rate uj is given by (12.11) , where Q is defined in (12.111) . The second term on the 
RHS of (12.21) is the rate at which cells transfer into A from Q and Q*. The transition rate a 
is given by (12. II) .and A is given by (12.111) . Only f2 and f2* cells in the G\ phase (i.e., with 
time counters c between c\ and C2) can transfer into A, which explains the boundaries in the 
integral in H 2 . 2 1) and the indicator function in (12.41) . In H 2 . 3 f) . the first term on the RHS is 
the rate at which cells flow from A into A* . These A cells flow from the endpoint x = into 
A*. The second term on the RHS of (I2.3P is the rate at which cells flow out of A* into Q*. 
Cells coming from A* enter Q* at the point (x, c) = (0, 0). 

The expression on the RHS of (\2A\\ represents cells in the Gi-phase that transfer out of 
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Q. The RHS of (12.51) accounts for the rate that cells in the Gi-phase flow out of Q* into A. 
Boundary condition (12. 6p indicates that all Omega cells that attain the maximal value of 
x = 1 commit to differentiation and do not transfer back to A. The first term of boundary 
condition (12. 7p indicates that once cells come to the end of a cell cycle at c = C2, they reset 
their time counters to c = 0. The second term accounts for the rate that cells transfer from 
A to Q at the boundary c = 0. It is the negative of the first term on the RHS of (12. 2p . 
Boundary condition (12. 8p accounts for the rate that cells transfer from A* to Q at the point 
x = 0. This expression balances the second term on the RHS of (12.31) . scaled by the advection 
rate because it represents a flux. Boundary conditions (12. 9p and (12.101) account for division 
whenever cells cross the point c = c\ of the cell cycle. 



2.1 Rescaled parameters 



For convenience, we rescale the parameters from the original ABM [35] to make them more 
appropriate for the PDE formulation presented in Section [2j (The values of the parameters 
used in |35j are listed in Table IA.1I ) First, we measure time in units of days rather than 
hours. Then, instead of using the affinity to characterize a stem cell's proclivity to remain 
in the Alpha (quiescent) state, we replace the affinity variable, a, with x, which is defined 
by 

a(x) = e 1X , 7 = -loga min = 6.2146. 

Since the affinity ranges from a m i n = 0.002 to a max = 1, x ranges from to 1. 
In [35], at every time step of 1 hour, the affinity increases by a factor r for cells in Alpha 
and decreases by a factor of d for cells in Omega. Hence, the rescaled advection rates are, 
taking into account the time unit of days, 

24 log r 24 log d 
p r = — = 0.3681, p d = — = 0.1884. 

7 7 

We also rescale the cell populations down by a factor of 10 5 cells, so that the scaling factors 
Na and Aq become equal to 1 rather than to 10 5 . Unlike |35j , we write the transition rates, 
a(A(t),x) and u(Q(t),x) as functions of x and not of a. Finally, since time is measured in 
days, we need to multiply the original transition rates f a /u by 24 as shown in equation (lA.ip . 
Our rescaled parameters are listed in Table IA~T1 alongside the original values from [35J. 



3 Approximations of the PDE system 



The PDE system of Section [2]is rather complex and hence difficult to analyze mathematically. 
Thus, we introduce several approximation steps that can be used to simplify the system. The 
approximation steps are as follows: 
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(0) No approximation. This label refers to the full PDE model (I2.2I) - (I2.10I) presented in 
Section [2J 

(1) Eliminate the explicit representation of the cell cycle by averaging the system over the 
variable c. 

(2) (a) Assume the transition rate u(x, Q) vanishes for all x, except x = 0. In other words, 
Alpha cells do not transition into Omega, except at the point (x,c) = (0,0). This 
approximation causes the variable Q vanish. 

(b) Assume the advection rate p r is effectively infinite so that cells transfering from 
Omega into Alpha immediately enter A* without passing through A. This approxima- 
tion causes the variable A to vanish. 

(3) Combine Approximations 1 and 2. 

(4) Combine Approximations 1 and 2 with the further approximation that the transition 
rates a(x, A) and Q(x, fl) do not depend on x. 



3.1 PDE system for Approximation 1 



In this approximation, we eliminate the explicit representation of the cell cycle by averaging 
the system over the variable c. which is removed, allowing us to combine Q and Q* into 
one single population. The new unknown Q now represents the total Omega population, it 
satisfies a new PDE that replaces equations (12.41) and (I2.5p . namely 

dQ dfl — — 

-5T + Pd^r = u{x, Q)A(x, t) + (-Ka(x, A) + b)Q(x, t). (3.1) 
at ox 

This averaging procedure leads to continuous rates of division and transfer, b and na(A,x), 
where b is the average growth rate and k is the proportion of time a cell spends in the 
Gi-phase of the cell cycle. These parameters are estimated in Section [331 
The boundary condition for Q at x = is the same as the original boundary condition for 
Q* given by (J2J|, i.e., 

n(0,t) = "Q'ty A*(t), where U(t) = f Q(x,t)dx. (3.2) 

Pd Jo 

We complete Equations (13.11) and ( 13. 2p with the following equations for A and A*. 
dA OA — 

— — P t ~q~ = ~ UJ ( x i Q)A(x, t) + Ka(x, A)Q(x, t), (3.3) 
dA* — 

— = Pr A(0,t)-u(0,tt)A*(t) 7 (3.4) 

A(l,t) = 0. (3.5) 
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To obtain the PDE (13. 3p for A we have just modified (12.21) to take into account the changes 
in the model for the Omega cells. The last term on the RHS balances the corresponding 
term in (13.11) . The ODE for A* and boundary condition for A remain unchanged (equations 
(13 .4p and (13.51) are exactly equations (I2.3P and (I2.6p . Approximation 1 is essentially the same 
PDE system used by Roeder et al. in [31]. 



3.2 PDE system for Approximation 2 



In Approximation 2, we suppose A and Vt are negligible compared to A* and Vt* . Indeed, 
numerical simulations from (23] show that over 98% of Alpha cells remain in the A* com- 
partment and over 94% of Omega cells remain in Q* compartment over time. Hence, to a 
good approximation, we can exclude the A and Q populations. As a result, we are only left 
with the following equations 



dt dx 
with boundary conditions 



-u{0,Q)A*(t) + f a{x,A*)l Gl (x)Q*(x,t)dx, (3.6) 
Jo 

f)C>* (jO* 

+ Pd — = -a(x,A*)l Gl (x)n*(x,t), (3.7) 



dA* 
~dt 



£T(0,t) = A*(t),tl(t) = [ ft*(x,t)dx, (3.1 



o 



n*( Pd (kc 2 + Cl ) + ,t) =2n*(p d (kc 2 + c 1 )-,t), k = 0,1,2,.... (3.9) 

The first term on the RHS of (I3.6P for A* accounts for the rate that cells flow from A* into 
Q*, and the second term on the RHS accounts for the rate that cells transfer from f2* directly 
into A* without passing through A. In the equations above, (13.61) replaces (I2.2p . (I2.3P and 
(12. 6p and can be obtained by integrating (12.21) with respect to x, adding it to (12.31) . and 
assuming that 

A* ^> / A(t,x)dx and / a(x, A)lo, 1 (x)Vt*{t,x)dx ^> I I a(x, A)Vt{t,x^c)dxdc. 

JO Jo Jo Ja 

Equation ( 13. 71) and boundary conditions f]3. 8[) and ( 13.91) for Q* remain the same as (12.5p . 
(EHD and (12~T0D . respectively. 
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3.3 PDE system for Approximation 3 



We combine Approximations 1 and 2 and to simplify notation, replace na(x, A) by a(x, A*) 
and u;(0, Q) by u(Q). Then we arrive at the system 

a A* _ r l 

— = -uj(n)A*(t) + / a(x,A*)n*(x,t)dx, (3.10) 
at J 

,9Q* 3Q* 

~dt + Pd ^x~ = A * ] + h) t} (3,11) 

with boundary condition (12.81) . i.e., 

n*(o,t) = ^@-A*(t), U{t)= [ n*(x,t)dx. (3.12) 

Pd Jo 

The parameter b is defined in the same way as in Section 13.11 and is estimated in Section 13.51 



3.4 PDE system for Approximation 4 



We combine Approximations 1 and 2 with the further approximation that the transition 
rate a(x,A) does not depend on x. To do so we may replace a(x,A) with its midpoint 
approximation a(A*) = a(0.5, A*). In addition, we write u>(f2) in place of cc;(0, Q). Thus, we 
arrive at the simplified system 

dA* _ r 1 

— = -u(Sl)A*(t) + a(A*) / n*(x,t)dx, (3.13) 
dt J 

,9Q* ,9Q* 

~dt + Pd ~dx~ = { ~ aiA * ] + b) t} ' (3 ' 14) 
with boundary condition 

fi*(0,f) = ^i*(t), U(t)= f Q*(x,t)dx. (3.15) 
Pd Jo 



3.5 Additional parameters 



The parameters b and k depend on how frequently cells circulate back and forth between the 
Alpha and Omega compartments. In particular, frequent turnover results in a decreased cell 
cycle. In |34j . Roeder et al. estimate b and k as functions of the total cell population, A + Q. 
For simplicity, we assume the total cell population remains constant at the value 1 x 10 5 . 
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Then using the functions in it follows that Omega cells divide once every 39.5 hours on 
average and that they spend 54% of their time in the Gi-phase |34j . Hence, we estimate 

b = = lQ S( 2 ) = o.42/days and k = ^— = 0.54 

ci + C2 39.5 hours c\ + C2 

These estimates are listed in Table IA.U 



4 Numerical simulations 



We compare numerical solutions of the four PDE models presented in Section[3]to simulations 
of the ABM from [31]. For the transport PDEs we use the explicit upwind scheme while the 
ABM uses the algorithm from [35J, which is summarized in Appendix [B] All programs are 
run in Matlab 7.6.0. 

Figure I4~T1 shows time plots of numerical solutions to Approximation (i.e., the full PDE 
model) and simulations of the ABM for 100 days. In each of the plots, the initial condition 
is A*(t — 0) = 1 while all other variables are 0. All parameters are taken from Table |A~T| 
except in Figures 14.1b . b, and c, where the differentiation rate d = 1.02, 1.05, and 1.2, 
respectively. (For the PDE model, these values of d yield p d values of 0.0765, 0.1884, and 
0.7041, respectively.) 




50 100 50 100 50 100 



Time (days) Time (days) Time (days) 

(a) (b) (c) 



Figure 4.1: Time plots of numerical solutions to Approximation and simulations of the 
ABM. Solutions to Approximation are shown by solid black lines, and simulations of the 
ABM are shown by dashed gray lines. The results of Approximation and the ABM are 
almost indistinguishable in all cases, except for the Omega population in plot (a), (a) Time 
plots for d = 1.02. (b) Time plots for d = 1.05. (c) Time plots for d = 1.2. 



The plots in Figures 14.1k . b, and c show the following characteristic behaviors, respec- 
tively: approach to periodic behavior, approach to a nonzero steady state, and approach to 
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the zero steady state. Furthermore, we see that the behavior of the full PDE model closely 
matches the behavior of the ABM in all three cases. 

Using the same initial conditions and parameters as above, Figure 14.21 compares Approx- 
imation 1 to the ABM. We see that the behavior of Approximation 1 qualitatively matches 
the behavior of the ABM to a large extent in all three cases. 




Figure 4.2: Time plots of numerical solutions to Approximation 1 and simulations of the 
ABM. Solutions to Approximation 1 are shown by solid black lines, and simulations of the 
ABM are shown by dashed gray lines, (a) Time plots for d = 1.02. (b) Time plots for 
d = 1.05. (c) Time plots for d = 1.2. 



Using the same initial conditions and parameters as before, Figure fl~3l compares Approx- 
imation 2 and the ABM. We see that the behavior of Approximation 2 closely replicates the 
behavior of the ABM in all three cases. In fact, the results of the two models are almost 
indistinguishable in all cases. 

Using the same initial conditions and parameters as before, Figure H~4l compares Approx- 
imation 3 and the ABM. Although both Approximations 1 and 2 followed the behavior of 
the ABM very well, in this case, we see that combining the two approximations leads to 
quite different behavior. 

Using the same initial conditions and parameters as before, Figure H~5l compares Approx- 
imation 4 and the ABM. Approximation 4 behaves a lot like Approximation 3, so the added 
assumption that the transfer function a does not depend on x does not affect the dynamics 
of the model very much. 
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50 100 50 100 50 100 



Time (days) Time (days) Time (days) 

(a) (b) (c) 

Figure 4.3: Time plots of numerical solutions to Approximation 2 and simulations of the 
ABM. Solutions to Approximation 2 are shown by solid black lines, and simulations of the 
ABM are shown by dashed gray lines, (a) Time plots for d = 1.02. (b) Time plots for 
d = 1.05. (c) Time plots for d = 1.2. 




Figure 4.4: Time plots of numerical solutions to Approximation 3 and simulations of the 
ABM. Solutions to Approximation 3 are shown by solid black lines, and simulations of the 
ABM are shown by dashed gray lines, (a) Time plots for d = 1.02. (b) Time plots for 
d = 1.05. (c) Time plots for d = 1.2. 
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Figure 4.5: Time plots of numerical solutions to Approximation 4 and simulations of the 
ABM. Solutions to Approximation 4 are shown by solid black lines, and simulations of the 
ABM are shown by dashed gray lines, (a) Time plots for d = 1.02. (b) Time plots for 
d = 1.05. (c) Time plots for d = 1.2. 



5 Analysis of Approximations 3 and 4 



In this section, we study both the system fl3.10l) - fl3.12l) corresponding to Approximation 
3 and the simpler system (I3.13p - fl3.15p corresponding to Approximation 4. To conduct 
the analysis we use duality arguments related to the "General Relative Entropy" method 
introduced in [31], [26] and used widely (see [TQl [91 [Tl [HI US] for other examples of applications). 
This method requires us to handle the eigenvalue problem and its adjoint, which we do first, 
and then use them to build entropy functionals. Another method is to reduce the system to 
a delay differential equation (DDE), and we also comment on how this is possible. 



5.1 Link with a Delay Differential Equation 



In the system (13. 130 - 03.1 51) . we can interpret Q* as the maturing and proliferating stem 
cell compartment, whereas A* represents a reservoir of quiescent, completely immature stem 
cells. There are exchanges between these compartments, with the same rules as before, 
except that once a cell stops the maturation process and enters the quiescent compartment, 
it immediately becomes fully immature, but if it reenters the maturing and proliferating 
compartment, it must go through the entire maturation process again. The main idea 
remains, however, unchanged, i.e., maturation is considered an invertible process as long as 
full maturity (at x — 1) is not reached. 

Although this assumption is original, and up to our knowledge has first been proposed 
in the ABM model of [22], the simplified system fl3.13p - fl3.15p . has some resemblance with 
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the CML models based on delay differential equations (DDEs), proposed by Mackey et al. 
[21EIHE3CE31E2]. Indeed, we can convert the PDE system (1333]) - (1535]) into an equivalent 
DDE system by applying the method of characteristics to (13.141) to obtain 

Q*{x, t) = Q* (o, t - — \ exp ( — x - [ a{A*{u))du 



Pd Jt-x/p d 
t 



A* (t- — \ exp( — x- [ a(A*(u))du\ 

V Pdj {Pd Jt-x/p A J 



Pd \ Pdj IPd Jt-x/p d 

for 1 < pdt, a restriction that discards the initial data and simplifies the setting. Integrating 
( 13.14p with respect to x and using the equation above to express Q(t, 1) and Q(t, 0) in terms 
of A*, we obtain 

^ + oo(U)e^~ C{t) A*(t - 1/ Pd ) - lu(U)A* = (b- a(A*)) O, 

where 

C(t) = J a(A*(u))du. 
Thus, we obtain the following DDE system: 

^ = u(t)A*(t) - u(t)e^- C{t) A*(t - 1/ Pd ) + (b- a(t))Q(t), (5.1) 
HA* — 

^— = -v(t)A*(t) + a(t)Q{t), (5.2) 

(Jjv 

dC 

— =a (t)- a (t-l/p d ), (5.3) 

where we use a(t) and u(t) to denote a(A*(t)) and u(Q(t)). The characteristic equation 
obtained from the DDE system ( 15. Ill — ( [5731 coincides with the eigenvalue equation ( 15.311) 
given by our analysis of local stability around the nonzero steady state of the PDE system 
(133T3|) - (I3~T5|) . (See Section EU) 



5.2 Existence of Steady States for Approximations 3 and 4 



The model based on Approximation 4 is simpler and we analyze it first. The method is 
extended to Approximation 3 afterwards. Steady states refer to time-independent solutions 
(A*, ft*) of system ( 13.13l) -( l3.15l) . They can be found and classified according to the model 
parameters. 
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Proposition 5.1 Let a(-) and u(-) be continuous positive bounded functions decreasing to 
zero. Then zero is always a steady state, and there is a nonzero steady state (which is unique) 
iff a(0) > b* , where b* is uniquely determined by b = b* for b = pd and by the relation 

b_ b* 

be Pd =b*e Pd : by^b*, for b ^ p d . (5.4) 
If a(0) < b* , then zero is the unique steady state. 



0.4 




1 1 1 1 1 1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

X 



Figure 5.1: The function F(x) = xe p d used in (15.41) . Its maximum is attained for x = pd- 

Proof. Obviously (Q = 0, A* = 0) is a steady state. We now consider possible nonzero 
steady states. Solving (13.141) . we find that nonzero steady state solutions (Q(x),A) satisfy 

b-a(A) 

Q(x) = n(0)e x , (5.5) 

where f2(0) ^ 0. Substituting (15.51) into H3 . 1 3|) and (13.151) . we determine the solution by the 
values (A, Q) that satisfy 

A = ^U, Pd = a{A) / e~^d~ x dx, Q = fl(0) / e~^T x dx. (5.6) 

Jo J 

Thus, we have two cases: 

Case 1. b = pa- This is equivalent to b = a(A) by the second equation of (15.61) . Therefore 
b* = b = a(A) and we can choose an A* = A ^ under the condition ct(0) > b* . Then the 
monotonicity condition on uj(-) allows us to find a unique Q by the first relation in ( 15. 6ft . 

15 



Case 2. b ^ p& and thus b ^ a(A). Since the second relation also implies 

_b_ _, a(A) 

be p d = a(A)e p d . 
we clearly have b* = a(A) and the result follows as before. □ 



Proposition 5.2 Let a(x, ■) and uj(-) be continuous positive bounded functions decreasing 
to zero (for all x in case of a). Then zero is always a steady state, and there is a nonzero 
steady state (which is unique) iff 

1 * b - a(yfi) „1 j cc(v,0)-b d 

a{x,0)e° Pd dx>p d ^ e* Pd dx>^. (5.7) 
I) Jo ° 



We leave to the reader that this condition is equivalent to a(y, 0) > b* in the case a is 
independent of stated in Proposition 15.11 

Proof. The explicit solution is now given by 

j b- a (y,A*) 

Therefore, as in the proof of Proposition 15. 1\ one reduces the claim to first finding A* that 
solves 

a ( X ' A *) c u(x,A*) d 
o Pd 

with u(x, A*) = J® b ~ a ^ A ' ] dy. Because a<yX p ^ - = -^—u' x (x, A*), integration by parts reduces 
our condition to 

1 = — [ e ^^) dx + e u{o,A') _ e u(i,A*) 



Pd Jo 

or also 



rl rl r a( y ,A*)-b , 

/ e <^)-u(i,A")d x = / e i '* dx=:H{A*). 
Jo Jo 



b 



This function H(-) is decreasing and satisfies 

H(oo) = ^(1 - e~^ b ) < ^. 



Therefore our condition (15.71) on H(0) is necessary and sufficient for the existence of a solution 
A* > 0. The end of the proof is the same as before. □ 
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5.3 An eigenvalue problem for Approximation 3 



Our analysis of a priori bounds and stability uses the eigenelements of the linear problem 
related to various states. This section introduces the relevant material. 

Here we fix A* and Q*, we set a(x) = a(x,A*) and uj = uj(Q), and we consider the 
eigenvalue problem to find (A G R, Q(x),A) that satisfy 

dfl 

Pd -^ = (b-\-a)n, 0<x<l, (5.8) 
0= / aVt(x)dx - (A + u)A, (5.9) 



o 



n(o) = —a, f 

Pd Jo 



Q(x)dx = 1. (5.10) 



Notice that we can simplify the problem by combining the last two relations to 

(X + uj)Q(0) = — [ att(x)dx, [ Q{x)dx = l. (5.11) 

Pd Jo Jo 



Proposition 5.3 Let a(x) > and uj > 0. There exists a unique solution X, Q(x) > 0, 
A > to A5.8\) - IR)7nfy . This X is the unique real eigenvalue, and it is positive iff 

1 J_ f(b~a(s))ds 

a(x)e" d o dx > p d , (5.12) 



and it is negative iff 

1 -i f(b~a(s))ds 

a(x)e" d o dx < p d . (5.13) 

o 

// a(x) = a is a constant, then A5.1fy) reduces to a > b* where b* is defined by \5.4\ )- The 
eigenvalue equals iff a = b*. 



Proof. One calculates explicitly from (I5.8P that 

-i- f(b-a(s)-X)ds 

n(x) = fi(0)e Pd o , 
and together with (15.111) . that the eigenvalue A must satisfy 

fX \ f 1 . . j- d h-a(s)-X)ds 

Pd - + 1 = / a(x)e dx. 



UJ 
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We can write this relation in the form G(X) = L(X), with 



G(X) = Pd [- + l 



L(X) 



1 / \ f(b-a(s)-X)ds 

a{x)e ax. 



The function G increases linearly from — oo to +00, whereas L(\) decreases continuously 
from +00 to 0. Hence, the two curves intersect at a unique value A G R. Moreover, A > 
iff G(0) < L(0), which gives Equation f)5.12p . If a(x) is constant, we obtain 

1 < [e fd - 1 

b — a \ 

which, with the notation F(z) = ze~ z ^ pd , simplifies to 

b — a b — a 

If b > a, this inequality is equivalent to F(a) > F{b), so a > min(6, b*) = b*. If b < a, the 
inequality is equivalent to F(a) < F(b), so a > b*. Thus, in both cases, a > b*. □ 

Remark. The adjoint problem of fl5.8l) - fl5.10l) is 



-Pd- 



(b - A - a(x)) <p + a(x)*, 



dx 

(A + cu)^ = cu4>(0), 
0(1) = 0, 

and its solution <ft can be explicitly calculated as follows: 



— L f(b-a(s)-X)ds 

(x) = 0(O)e Pd o 



1 _ Jo 



a(s)e ds 



a(s)e Pd ° ds) 



(5.14) 

(5.15) 
(5.16) 



(5.17) 



5.4 Uniform Bounds 



The aim of this section is to show that the population remains bounded under the previous 
assumption that a vanishes at infinity (this can be relaxed as shown below). This relies on 
the eigenvalue problem and the property that the eigenvalue is negative for large populations. 



Proposition 5.4 If cj(-) is bounded, a(A,x) is a nonin creasing function of A and for A 
large enough 



, : , b- a (y,A) d 

J Pd ^ 



/ a(x,A)e° Pd "dx<p d , (5.18) 
Jo 
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then any solution (Q*(t, x), A*(t)) of nUKl - HJE) remains bounded for all t > 0, i.e., 
A* E L°°(0, oo) and ft* G L°°((0, oo) x (0, 1)). 



Proof. Since ( 15.181) holds, we can choose A large enough so that the condition (I5.13P is 
fulfilled for a(x) = a(x, A). According to Proposition 15.31 there exists a unique eigenvalue 
A < and a unique positive couple (A, ft) that solves fl5.8l) -( l5TT0l) with a(x) = a(x, A). We 
also denote by the solution of the adjoint eigenproblem f l5.14l) - fl5.16l) . Testing the solutions 
Q*(t,x) of Problem fl3TT0|) - fl3TT2|) against gives 



A 



J (j)VL*dx = J(f)(b- a(x, A*))n*dx + J p d n*£(f)dx + p d Sl*(t, 0)0(0) 

= fQ*(Xcf>-a(x,A)^)dx + f(a(x*A) - a(x,A*(t)))Q*(f>dx 
+p d ft*(i, 0)0(0) 

= £ ^3f i ^°)) + {a(x,A)-a(x,A*(t)))n*<i S jdx 

+0(0) (-TT + J? A * (m*dx) 
= /ft*(A0 + X^0^)dx + j(a{x,A) -a{x,A*{t)))Q*<j)dx 

+0(0) (-^ + / (a(x, A*(t)) - a(x, A))tt*dx) 

= JQ*(X4> + X^^)dx + j(a{x, A) -a(x,A*(t)))n*((j)(x)- (f)(0)) dx 

-mi* 

Hence, we have obtained the following equality: 

+ In* (0(x) - 0(0)) (a(x, A) - a{x, A*(t))) dx. 



We consider the quantity S(t) = f (p(x)Q*(x,t)dx + 4>(Q)A*(t) - pmm yA,A*(t)j , where \i 
is a constant such that 

Vxe[0,l], pa(x,A) > (-A + l)0+||a(.,.)|| oo ||0|| oo . (5.19) 

Recalling that A < 0, we obtain 

^— < X I ft*0 + I ft* (cj)(x) - 0(0)) (a(x, A) - a(x, A*(t)fj dx 

+ pl A * mA (uA*(t) - J a(x,A*( y t))Q*( y t,x)dx\ . 

We set 5*o = 0(O)^4 + /i||u;|| OO 74 and we prove that S{t) < max(S'(0), So). Indeed, If S(t) > So, 
we have two cases: 
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First case: A*(t) > A. The last term in the RHS of (15.201) vanishes, the first two terms 

dS 
dt 



are less than or equal to zero, so < and S(t) decreases. 



Second case: A*(t) < A. Equation (15.201) implies (recall that a(x, A*) > a(x, A) ) 
'' S <\ I ft*4>(0)a(x,A*(t)) +fjwA*(t) -/I [ a(x,A*(t))Q*(t,x)dx. 



dt 

We deduce, using (15.191) . that 

dS 



o 



, < fiWuWooA - I n*(f)dx < S - S(t) < 0. 

(JjV 







This shows that S(t) decreases in both cases, which proves that S(t) < max(5*(0), Sq) 
remains constantly bounded. As a consequence J <ft(x)Q*(x,t)dx + <f)(0)A*(t) is uniformly 
bounded. Then, since A*(t) is bounded, in (I3.12p the boundary term f2*(0, t) is also uniformly 
bounded and thus Q* e L°°((0, oo) x (0, 1)). n 



5.5 Stability Analysis near Zero 



Proposition 5.5 Ifa(x,A) is a positive decreasing function of A, and if 

1 -A- f(b-a(s,0))ds 

a(x,0)e Pd « dx< Pd , (5.21) 



o 







then the zero steady state is globally attractive. On the contrary, if 

a(x, 0)e^ fo(»-^mds dx > p ^ (5 22) 

then the zero steady state is unstable. 



Remark. According to the proof of Proposition 15. 2\ the stability condition (15.221) is equiva- 
lent to condition (I5.7P for the existence of a nonzero steady state. Therefore the zero steady 
state is stable exactly when it is the unique steady state. Also, the properties of u;(0, Q) do 
not influence the stability of the zero steady state. 

Proof. We consider the linear eigenvalue problem (I5.8I) - (I5.10I) at zero and let 0, given by 
(I5.17p . be the solution to the adjoint eigenvalue problem. Let (fl(t, x), A(t)) denote a solution 
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to the time-dependent problem (I3.13I) - (I3.15I) . As in the proof of Proposition I5.4[ we test Q 
against and define v(t) = J <p(x)Q(x,t)dx + 4>(0)A(t). We have 

dt J V ) (5.23) 

+ / n (0(x) - 0(0)) (a(x, 0) - a(x, A)) dx. 



If A < 0, which is equivalent to (15.211) according to Proposition 15.31 (15.171) shows that is 
decreasing, and since we have supposed that a decreases with respect to A, we can conclude 
that 

d 



dt 



4>Qdx + (f>(0)A < A J <pttdx. 



It follows that J <j)Qdx + 4>(0)A is decreasing. Since this integral is nonnegative, it tends to 
a limit, and at infinity, f <pQdx tends to 0. Therefore, from the uniform bounds, Q(t) — > 0, 

and thus A(t) — ► 0. Thus, the zero steady state is globally attractive. 



t— >oo 



If A = 0, we still have that v(t) is nonincreasing, and one has 

—v(t) = / Q(x, t) (0(x) - 0(0)) (a(x, 0) - a{x, A))dx. 

Since <p(x) — 0(0) < and a(x, 0) — a(x, A) > for x < 1 and A > 0, it implies that either 
Q(x < 1, +oo) = or A(+oo) = 0. We conclude using fl3TT3D -f l3TT5l) that Q = A = 0. 

If A > 0, we define the constant C = ||0'||oosup \a' A (x,A)\/</>(0) where the sup is taken 
on < x < 1 and < <f){0)A < oo. Then, we write 



d r l 

—v{t) > / 0fi[A - C<f>{0)A(t)]dx 
dt J 



> / (pQ[X-Cv(t)}dx. 
Jo 

This proves that v(t) is increasing whenever v(t) < X/C and the zero steady state is thus 
unstable. □ 

Remark. We have not made any assumption on Indeed, the population in the A* 

compartment neither divides nor dies. Hence, whether cells in A* transfer quickly or slowly 
into the Q compartment does not change the behavior at infinity or at zero. On the contrary, 
we can expect it influences the stability of the nonzero steady state, whenever it exists. 
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5.6 Stability analysis near the nonzero steady state 



To study the local stability of the operator around the nontrivial steady state, denoted Q, A, 
we consider a small perturbation 5Q = Q — Q, 5 A = A — A, where A(t,x) and Q(t,x) are 
solutions of System H3. 13[) — H3. 15[) . 

For the sake of simplicity, we denote f2 = j Qdx, a = a (A) , a' — ^ (A) , cu = cu (fi) and 



UJ 



I dw 



At first order in 5 A and 8Q, the perturbation 5Q, 5 A satisfies the following equation: 

d5Q dSQ . ,~ . , nAS 

~df Pd ~dx~ = ( b -^ 5n - a QSA i ( 5 - 24 ) 

—- = («- J A) / *n(*, x)cfe + ( / a'Q(x)dx - u)5A, (5.25) 
dt Jo Jo 

8Q(0,t) = —6A(t) + — [ 5Q(t,x)dx. (5.26) 
Pd pd Jo 

We now consider the eigenvalue problem related to the linear system f l5.24p - fl5.26p . Using 

b-b* 

ct(A) = b* = be p d (cf. Equation (JH3J) and the definition of Vt in f)5.5p . we calculate 

n(x)dx = n(o) Pd = — — . 

We also recall that Q(0)pd = Auj by (15.61) . Hence, the eigenvalue problem writes 

A5ft + Pd — = (b- b*) 6n - c/&Afi(0)e~^~*, (5.27) 
ax 

AM = b* f 5Q(x)dx + a'5A^P± _ ^K^^ _ ^ /" SQ(x)dx, (5.28) 
Jo "* ^4 Jo 

SQ(Q) = ^5A + — C 5Q(x)dx. (5.29) 



A pd Jo 

(We still denote the eigenvector by (5 A, 5Q).) By explicitly solving f)5.27p . we obtain 

Sn(x) = e~^ x I 5n(0) - a'tt(0)5A- 



X 



We impose the restriction 



/ 5Q(x)dx = 1 
Jo 



<K](0)p d , 6 _a a'O(0)^U(l-e^)-A 

7T e Pd - !) \ II E £ x — , (5-30) 



b-b* - \ x b* ' A b* b-b* - A 
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and obtain from f 1 5 . 2 8 1) that 



X5A = b* + a '5A^- Pd - ^-p d 5A - J A. 

o* A 



In addition, (I5.29f) yields 



H A X- a ^p d + m Pi 
which we substitute into (15.301) to obtain 

1 = /(A):= 




uj'A 



(5.31) 

a^(0)p d o ^- e ~ A b*-m'A 

b-b*-\ x _ a ,m Pd + m Pd 

Relation (15.311) defines all eigenvalues of the linearized system (I5.27I) - (I5.29I) . It is the same as 
the characteristic equation obtained from the DDE system (I5.ip - (I5.3I) . The transition from 
stability to instability regions has been studied in many cases of delay equations and we can 
expect that it corresponds to a Hopf bifurcation that explains the appearance of periodic 
solutions, see [151 [25]. Section [6] numerically studies this equation, which shows that the 
nontrivial steady state can either be stable or unstable depending on the parameters of the 
equation. 



6 Stability regions 



Since the analytic expression of the eigenvalues in Section 15.61 is difficult to work out, we 
determine them numerically. To do so we use the Matlab program DDE-BIFTOOL[18] to 
determine, equivalently, the eigenvalues of the DDE system (I5.1l) - (l5.3p as the parameters pd 
and b vary. In this way, we can determine the regions of stability for the nonzero steady 
state. Futhermore, applying Proposition 15.51 we can numerically determine where this state 
disappers and the zero steady state becomes stable. 

Notice that this procedure is equivalent to using the eigenvalue equation (I5.3ip to nu- 
merically determine stability regions for the PDE system (I3.13p - (l3.15p . However, since the 
DDE-BIFTOOL software already exists, it is easier to numerically analyze the equivalent 
DDE system given by (I5.1l) - (l5.3p instead. 
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Figure [67T1 shows the stability regions of the system with respect to pd and h. The variable 
pd denotes the advection rate of Omega cells with respect to x, measured in units of 1/day. 
The variable b is the exponential growth rate of Omega cells, measured in units of 1/day. 



7r 




0.2 0.4 0.6 0.8 1 



Pd 

Figure 6.1: Stability regions for the PDE system from Section [3.41 in (pd, 6)-space. The 
variable pd denotes the advection rate of Omega cells with respect to x, measured in units 
of 1/day. Higher pd means that cells spend less time in the Omega compartment before 
differentiating. The variable b is the exponential growth rate of Omega cells, measured in 
units of 1/day. (a) In the top region, the system is unstable and exhibits periodic behavior, 
(b) In the middle region, the system is stable at the nonzero steady state. Hence, the Alpha 
and Omega populations approach a nontrivial equilibrium, (c) In the bottom region, the 
nonzero steady state does not exist, and the system is stable at the zero steady state. 

As we can see from Figure I6.1[ the region at which the system is stable at the nontrivial 
steady state falls between the regions where the system is unstable and where the system is 
only stable at zero. Hence, the stability of the system at a nonzero equilibrium depends on a 
balance between the rate of differentiation and the growth rate of the Omega (proliferating) 
stem cell population. 

If the stem cell growth rate increases or if the rate of differentiation decreases sufficiently, 
the system transitions to instability, in which case population exhibits periodic behavior. 
This transition from a stable equilibrium to unstable periodic behavior might correspond to 
the transition between CML and acute myelogenous leukemia (AML), a disease characterized 
by the rapid proliferation and invasion of the blood by immature, undifferentiated cells. 
Indeed, unstable and high amplitude oscillations in the stem cell population would result in 
the sudden overproduction of immature, progenitor cells. These results concur with those 
of PP, which presents a DDE model of hematopoiesis and concludes that obstructing cell 
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different at ion at an early stage of development results in the overproduction of immature 
cells, potentially leading to AML. 

On the other side, if the growth rate decreases or the differentiation rate increases suf- 
ficiently, the system loses the nontrivial equilibrium and becomes stable only at zero. This 
result makes sense, since the stem cell population has to proliferate fast enough to replenish 
itself as Omega cells continually differentiate. However, it is currently unclear to us whether 
this scenario corresponds to a particular disease. 

Finally, the program DDE-BIFTOOL can also be used to numerically determine the 
trajectories of the rightmost eigenvalue of (I5.3ip as it crosses the imaginary axis from stability 
to instability (i.e., from negative to positive real part) as parameter pd or b is varied. The 
trajectories of the rightmost eigenvalues as pd and b are increased independently are shown 
in Figure I6.2L As also shown in Figure 16.11 we see that the system moves from instability to 
stability as pd increases, whereas the system moves from stability to instability as b increases. 




i , , , , 1 , i , 1 ^ 

-2.5 -2 -1.5 -1 -0.5 0.5 -1 -0.5 0.5 



(a) (b) 



Figure 6.2: Trajectories of the rightmost eigenvalue around the first stability crossing point. 
Each trajectory is associated with a conjugate trajectory (not shown), (a) Trajectory of the 
rightmost eigenvalue as pd varies from 0.0422 to 0.3505. The value of b is fixed at 0.42. (b) 
Trajectory of the rightmost eigenvalue as b varies from 0.2 to 1.5. The value of pd is fixed at 
0.1884 (d = 1.05), which is the estimated value from [35] 



7 Conclusion 



In this paper, we simplify the PDE model of hematopoiesis presented in [23] and perform 
a stability analysis. The PDE model from [23] is a time-continuous extension of the ABM 
of [35] and has shown to coincide closely with the dynamics of the original ABM. Our 
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simplifications render the original PDE model amenable to analysis while closely preserving 
the qualitative behavior. 



Prior to our analysis, we simplified the PDE model given by fl2.2p - fl2.5l) and fl2.6l) - fl2.10l) 
to the system labelled Approximation 4 in Section 13.41 As shown in Figure 14. 5[ Approx- 
imation 4 is a good approximation of the original ABM from a qualitative perspective. 
Furthermore, Approximation 4 can be analyzed analytically. Hence, the process of simplifi- 
cation and analysis in this paper, allows us to apply an analytical approach to an otherwise 
inherently complex and intractable ABM. Furthermore, the observation that Approximation 
4 is equivalent to a DDE system shows that the ABM of [35] possesses inherent similarities 
to other age-structured DDE models proposed by Mackey et al. [2], [3j HI [T2J, [131 E2] and 
Adimy et al. [I]. 



In addition, Approximation 2, presented in Section 13721 generates dynamics that are sur- 
prisingly close, even quantitatively, to those of the original ABM. (See Figure fl~3l ) Although 
Approximation 2 is not as suitable for analysis as Approximation 4, it possesses the same 
computational complexity, since it is also a system of two PDEs with one state variable. 
Indeed, it only takes minutes to run numerical simulations of Approximation 2 using the 
explicit upwind scheme. This speed of processing is comparable to that of Approximation 
4 and of the difference equation model presented in [22] ■ Hence, for the purposes of nu- 
merical simulations, Approximation 2 could be used as an alternative to Approximation 
4. In addition, since Approximation 2 is a PDE system, one could readily add additional 
ODEs, DDEs, or PDEs to the existing system to capture the concurrent dynamics of the 
anti-leukemia immune response during imatinib treatment as was done in [21] for the ODE 
model of imatinib dynamics in [28J. Since the dynamics of leukemia cells differs between the 
ABM model of [35] and the ODE model of [28], especially in the long-term time scale, it will 
be useful to examine the impact of the anti-leukemia immune response on the ABM as well. 
We leave this for a future work. 

Finally, from our analysis, we observe that the system exhibits three types of behavior: 

(a) Periodic behavior, (b) A stable nonzero steady state, and (c) Stability at zero. State (c) 
corresponds to a state in which the stem cell population cannot maintain itself, and hence 
dies out to zero. It would correspond to a condition characterized by insufficient blood 
production. State (b) corresponds to the desirable state, in which the hematopoietic stem 
cells remain at equilbrium and can maintain a stable population level of blood cells. State (a) 
corresponds to a state in which stem cell and blood cell populations fluctuate rapidly. While 
state (b) captures the behavior of more stable blood cell growth such as in the case of CML, 
state (a) captures some of the relevant dynamics of AML. Hence, the transition from state 

(b) to state (a) might provide insight into the specific malfunction (either in differentiation 
or growth) that leads to the transition from CML to AML. 
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A Parameter estimates 



The sigmoidal transition functions are given in [35] by 



f a/u (A/0) = 24 ( + u A ] , (A.i; 

z/i + z/ 2 exp ( ^t- 7 "" 



where A and Q denote the total populations in the Alpha and Omega compartments, re- 
spectively (see (12. lip . Furthermore, 

hih 3 -hl fhs-vA 

u i = I I u oT - ' ^ = hi - Ui, u 3 = In , i/ 4 = / a / w (oo), 

h + h 3 - 2h 2 V ^2 / 



/a/«(0)-/«/ w (oo)' 2 f a / u (N A/{1 /2) -/ a/w (oo)' 3 f a / u (N A/a ) -/ a / w (oo)" 
The values of the various parameters are listed in Table IA.1I 



B Algorithm for the agent-based model 



We summarize the algorithm of the ABM from [35] . At every time step (1 hour), the ABM 
is defined as the following set of actions: 

A. Preliminary calculations. 
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Par am 


Description 


Ph" 


Ph + /imatinib- affected 


Rescaled 




Min value of affinity a 


0.002 


0.002 




flmax 


Max value of affinity a 


1.0 


f.O 




d 


Differentiation coefficient 


1.05 


1.05 


Pd = 0.1884 


r 


Regeneration coefficient 


l.f 


l.f 


p r = 0.3681 


Cl 


Duration of S/G2/M-phases 


17 hours 


17 hours 


17/24 




Cell cycle duration 


49 hours 


49 hours 


49/24 


A p 


Lifespan of precursor cells 


20 days 


20 days 


20 


A m 


Lifespan of mature cells 


8 days 


8 days 


8 


T c 


Division period tor precursors 


24 hours 


24 hours 


1 


fa(0) 


Transition characteristic for f a 


0.5 


f.O 




f*(N A /2) 


Transition characteristic for f a 


0.45 


0.9 




UN a) 


Transition characteristic for f a 


0.05 


0.058 




/q(oo) 


Transition characteristic for f a 


0.0 


0.0 






Scaling factor 


10 5 


fO 5 


1 




Transition characteristic for f u 


0.5 


f.O / 0.0500 




UN A /2) 


Transition characteristic for f u 


0.3 


0.99 / 0.0499 




UN A ) 


Transition characteristic for / w 


O.f 


0.98 / 0.0498 




U(oo) 


Transition characteristic for f u 


0.0 


0.96 / 0.0496 




Nn 


Scaling factor 


10 5 


10 5 


1 


Hnh 


Inhibition intensity 




0.050 


1.2 


r dcg 


Degradation intensity 




0.033 


0.8 


6 


Expansion rate for Approximations 






0.42 




Fraction of time spent in Gi-phase 






0.54 



Table A.l: Parameters from [35J and rescaled values used in this paper. The inhibition 
intensity, r in h, refers to the probability that a proliferative Ph + cell (i.e., an Q or Q* cell) 
becomes imatinib-affected in a given time interval. The degradation intensity, r^g, refers to 
the probability that an imatinib-affected, proliferative Ph + cell dies in a given interval. In 
the column for rescaled parameters, if the entry is blank, the corresponding parameters are 
left unchanged. 
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1. Calculate the total populations of A and Q cells. 

2. During imatinib treatment: 

• Remove the proliferative Ph+ cells and that undergo apoptosis. 

• Determine which unaffected proliferative Ph+, Q + , become imatinib- affected. 

B. Proliferation, death, change of state, clocks. At this stage, all cells fall into one of 
three categories: A stem cells, f2 stem cells, differentiated cells. 

1. For each A stem cell: 

• Determine whether the cell transfers to Q. If a cell transfers, skip the remaining 
actions for A cells. Note that the transition function depends on whether the cell 
is Ph~, Ph + , or imatinib-affected. Calculate transition probabilities based on the 
total population of f2 calculated in Step Al. 

• Increase the cell's affinity by a factor of r. 

2. For each Q stem cell: 

• Determine whether the cell transfers to A. If a cell transfers, skip the remaining 
actions for Q cells. Calculate transition probabilities based on the total population 
of A calculated in Step Al. 

• If the cell's affinity is less than or equal to a m i n , the cell becomes a differentiated 
cell of age 0. If the cell differentiates, skip the remaining actions for Q cells. 

• If a cell's affinity is greater than a m i n , decrease the cell's affinity by a factor of d. 

• Increase the counter c by 1. 

• If the counter c is greater than or equal to 49, set c to and create a new cell 
with identical attributes and state values as the current cell. 

3. For each differentiated cell: 

• Increase the cell's age by one. 

• If the cell's age is a multiple of 24 between 24 and 480, inclusively, create a new 
differentiated cell with with the same age as the current cell. 

• If a cells age reaches 672, that cell dies. 

Note that differentiated cells of age less than 480 are considered to be proliferating precur- 
sors, whereas differentiated cells of age greater than or equal to 480 are considered to be 
nonproliferating mature cells. 
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